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Abstract 

An accurate equation of state of the one component plasma is obtained in the low coupling regime 
< r < 1. The accuracy results from a smooth combination of the well-known hypernetted chain 
integral equation, Monte Carlo simulations and asymptotic analytical expressions of the excess 
internal energy u. In particular, special attention has been brought to describe and take advantage 
of finite size effects on Monte Carlo results to get the thermodynamic limit of u. This combined 
approach reproduces very accurately the different plasma correlation regimes encountered in this 
range of values of T. This paper extends to low T's an earlier Monte Carlo simulation study devoted 
to strongly coupled systems for 1 < T < 190 (J.-M. Caillol, J. Chem. Phys. Ill, 6538 (1999)). 
Analytical fits of u(T) in the range < T < 1 are provided with a precision that we claim to be not 
smaller than p = 10~ 5 . HNC equation and exact asymptotic expressions are shown to give reliable 
results for u(T) only in narrow T intervals, i.e. < T < 0.5 and < T < 0.3 respectively. 
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I. INTRODUCTION 



The aim of this paper is to obtain the equation of state (EOS) of a plasma in the low 
coupling regime with a high precision. In this regime standard Monte Carlo (MC) and 
Molecular Dynamics simulations techniques must be handled with care due to huge finite 
size effects and, in the other hand, the ideal gas approximation or more elaborated analytical 
expressions commonly used are valid only but asymptotically, for very small values of the 
coupling parameters. Such thermodynamic conditions are relevant for many astrophysical 
or laboratory plasmas hydrodynamics applications. 

However we shall restrict ourselves to the well known one-component plasma (OCP) 
model, which consists of identical point ions with number density n, charges Ze, moving in 
a neutralizing background, electrons for instance, where n = N/Q, N number of particles, Q 
volume of the system[l|]. In the very low coupling regime, the virial expansion supplemented 
by well documented resummation methods, as the well-known Debye-Huckel (DH) theory 
and its extensions (see e.g. Cohen 3J and, more recently, Ortner 4J expansions for instance) 
give reliable results. In the low to intermediate coupling regimes the HyperNetted Chain 
(HNC) integral equation ^] must be solved numerically. Finally, in the strong correlation 
regime, the OCP has also been extensively studied by Monte Carlo and Molecular Dynamics 



simulations for three decades, see e.g. 



, l8|, [9|, lid . LL3J and references cited herein. 



In the more recent of these references one of us has determined the thermodynamic 
limit of the excess internal energy per particle un=oo of the OCP with a high precision 
by means of MC simulations in the canonical ensemble within hyper spherical boundary 



conditions 



10 . [ll ] for 1 < T < 190. We recall that in the thermodynamic limit, i.e. for an 



infinite system of particles, the thermodynamics properties of the model depend solely on 
the coupling parameter T = f3(Ze) 2 /di (j3 = 1/kT, k Boltzmann constant, T temperature, 
and <2j the ionic radius defined by 47ma?/3 = 1), whereas, for a finite sample, an additional 
dependance on the number of particles N remains. In paper henceforth to be referred 
to as "I", special attention has been brought to describe and take advantage of such finite 
size effects on the energy un(T) to get its thermodynamic limit ujy =OQ , using all facilities of 
work stations available at that time. 

Recently we have also performed extensive MC simulations of the related Yukawa One- 
Component Plasma (YOCP), i.e. a system made of N identical point charges Ze interacting 



via an effective Yukawa pair-potential v a (r) = (Ze) 2 exp(— ar)/r, where a is the so-called 



screening parameter 



121 ]. not to be discussed however in this work. For the OCP and the 



YOCP as well, in the low T regime, the Debye length (i.e. the correlation length associated 
with charge fluctuations) becomes of the order or much larger than the size of the simulation 
box, yielding huge finite size effects on un(T). Therefore, despite these numerous studies 
and amount of work it appears that hydrocode applications using the combination of data 
bases and fits coming from various techniques can be affected by numerical instabilities in 
the transition regime, around T = 1. With nowadays computers it is now possible to explore 
this range of small T values with the help of performant simulation techniques and to obtain 
such precise results so that they can be considered as the reference ones to be used in many 
applications dealing with degenerate astrophysical or laboratory plasmas. We also examine 
carefully in this paper the connection between MC and first principle analytical or HNC 
results for T < 1. We have thus explored and precised the domain of validity of each of 
these methods. It turns out to be necessary to combine all of these approaches to obtain 
a continuous representation of un =00 (T) in the range < T < 1. Finally we extract from 
these combined approaches the best possible analytical representation for Uoo(r). 

Our paper is organized as follows. Next section is devoted to a brief presentation of 
the main features of low T expansions (Section III AO . the HNC integral equation (Section 
III Bj) and the rather unusual but efficient MC technique used in this paper (Section III CI) . 
Note that we have redone, by passing, extremely accurate HNC calculations and obtained 
new fits of HNC data, presented in Section III Bl In Section II II I we present and discuss our 
MC simulations. Fits of the data are described in details and widely illustrated. Finally 
conclusions are drawn in Section IIVI 



II. LOW r CALCULATION METHODS 



The interval < F < 1 covers various correlation regimes from no correlation (r = 0, 
i.e. the ideal gas) to an intermediate correlated regime (r = 1, no oscillation or structure 
in the pair correlation functions). In any case, the long-range nature of the interaction 
potential between two ionic charges causes Mayer graphs to diverge A field theoretical 
diagrammatic representation of cluster integrals has been proposed recently in [4] to avoid 
complicated chain resummations in an attempt to treat the T expansion of the classical 
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FIG. 1: Reduced excess energy ftujT versus < V < 1. Diamonds: HNC, black solid line: DH, 
thick cyan solid line: Thl approximation (|2.1aj) . thick cyan dashed line: Th2 approximation (|2.1b|) . 
other curves represent the successive orders of expansion (12. lj) . 



Coulomb system in a more controlled and systematic way. In this interesting paper the final 
expansion obtained by the author improves earlier and seminal analytical results of Cohen et 
al. 2j, y) obtained by traditional diagrammatic expansions and resummations. From these 
theoretical analysis it turns out that the physics in this small interval < T < 1 is extremely 
complicated and exhibits many different correlation regimes, even more than in the widely 



studied region 1 < T < 190 



10 



The low T expansions obtained by Cohen 



et al. and Ortner for Moo(r) converge to the HNC results only for < T < 0.2 as apparent 
in figure HJ For higher values of T these asymptotic expressions do not seem to converge 
at all and, moreover, the high order terms of the expansions do not improve the results of 
the lower orders. Anticipating the results of sections III Bl and III CI and, as can be observed 
in figure El the HNC data deviate from our MC results as soon as T > 0.5. It results from 
this sketchy discussion that we must distinguish three different regimes of correlations in 
the interval < T < 1, and we confess that this complexity motivated the present study. 
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FIG. 2: Ratio of MC excess energies to HNC results versus T in the low coupling regime. 
A. Cohen and Ortner analytical expansions 

In ref. {4] Ortner has developed an effective method based on the Hubbard- Stratonovich 
(HS) transformation and field theoretical approaches to calculate the free energy of classical 



Coulomb systems in the low F regime 13J, ll4j, Il5| . The HS transform was used to obtain 



the EOS of a classical plasma and notably that of the OCP. The non-trivial part of the 
Helmholtz free-energy density was derived up to order T 6 , improving on the previous results 

g 

of Cohen et al. at order Ta, obtained by a method of resummation of diverging diagrams. 
The author gives an analytical representation of the excess internal energy (3u of the OCP, 
valid at low T, without however any estimation of the error. It reads as, 



(3u{T) = p T 3/2 + Pl T 3 In T + p 2 T 3 + p 3 T 9/2 In T + p A T 9/2 

+ P5 T 6 in 2 r + p 6 r 6 in r + P7 v 6 



(2.1a) 
(2.1b) 



with the constants, p = —y/3/2, p\ = —9/8, p 2 = — (91n3)/8 — 3(7^/2 + 1, p 3 = 
-(27 v / 3)/16, p A = 0.2350, p 5 = -81/16, p 6 = -2.0959, p 7 = 0.0676 and C E = 0.57721566 
the Euler constant. Expression 12.11 (to be referred to as Th2 henceforth) improves on that 
given by Cohen et al. (to be referred to as Thl henceforth) [3], which corresponds to line 
I2.1al while the additional terms are those of line 12. lbl We recognize that the first term 
(— y / 3r 3 / 2 /2) is exactly the well known Debye-Huckel (DH) contribution. Figure [1] dis- 
plays the results of the reduced excess energy (3u/T versus Y at successive orders in the 
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FIG. 3: Reduced excess energy (3u/T versus V. Squares: MC (the symbols are larger than error 
bars), black line: DH, red line: HNC, blue line: Thl approximation (|2.1aj) . green line: Th2 
approximation (|2.1bj) , 

T-expansion 12.11 A close examination of the figures reveals that the DH approximation is 
nearly exact up to V = 0.05, in the sense that higher order contributions do not change the 
result. A comparison with HNC results, which are supposed to be nearly exact at least up 
to T = 0.5 (this point will be fully discussed in next section), shows the convergence of the 
expansions Thl and Th2 to HNC at T < 0.3 and T < 0.2 respectively. However we do not 
observe any trend of convergence of these expansions for T > 0.4. We also notice that the 
additional terms given by Ortner (cf equation I2.1b|) lead to an oscillatory behavior rather 
than to an improved convergence radius. We suspect some misprints in the reported p n for 
n = 5, 6, 7 since the T functional T dependence of 12.11 is undoubtedly correct. 

B. HNC method and fits 

1. Method 

We have redone high precision HNC calculations for a hundred of values of V in the range 
(0, 1) (see figures [H and [3] ) ; additional calculations were also done for some higher values of 
the coupling parameter, in the range 1 < T < 10, see figure HI We used the Ng method^ 
with the following control parameters: the pair correlation functions (direct and non-direct 
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FIG. 4: Comparison between HNC (red line) and MC data (squares, present work and previous 



results, see ref. 



lll |) for the reduced excess energy fiU/T versus 1 < Y < 10. 



respectively) c(r) and h(r), as well as their Fourier transforms c(k) and h(k), were tabulated 
on grids of iV = 2 M points with M = 20 in order to make use of fast Fourier transforms with 
intervals of Ar = 0.001 and Ak = 2tt/N ~ 610~ 3 in direct and Fourier space respectively. 
The dimensionless energies were computed according the formulae [l| 

pu^ 3 ro ° 



r 



dr rh(r) , 



dk h(k) 



(2.2a) 
(2.2b) 



r 2 (2tt) 2 J0 

where the distances "r" are measured in the units of the ionic radius cij and the wave numbers 
k in units of a" 1 . The comparison of these two estimations u^ r ' and of the energy, which 
of course should be equal, give an idea on the relative precision of the numerical resolution 
of HNC, typically about 10~ 12 at T = 0.01 and 10~ 13 at T > 0.1. Another usefull test is to 
check the Stillinger-Lovett (SL) sum rules ;recall briefly the two first SL rules (the third one 
should not be satisfied by HNc[l|) 



3 / dr r 2 h(r) 
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dr r 4 h(r) 
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r 



(2.3a) 
(2.3b) 



With the control parameters given above the SL rules were satisfied with a relative precision 
of about 10" 13 . 



TABLE I: First five Cohen-Ortner coefficients (cf Eq. (|2.ip . first line) compared to the correspon- 
dent coefficients of the fits of the energy (5u{T)/T for HNC and MC data. Second line : HNC, 7 
parameters, pO = — V3/2 fixed to its DH value. Third line : HNC, 8 parameters. Last line : MC 
data in the range 0.4 < T < 1 , 5 parameters (ps = pe = p-j = 0). 



pO 


pi 


p2 


P 3 


p4 


Method 


-0.8660254038 


-1.1250000000 


-1.1017662315 


-2.9228357378 


0.2350000000 


Ortner 


-0.8660254038 


-1.1127645260 


-1.0636075255 


-3.1960177420 


-1.4236810385 


HNC - DH 


-0.8658509448 


-1.0967358264 


-1.0224523661 


-2.9765709164 


-1.1861133643 


HNC 


-0.8409025523 


-0.5198391670 


-0.0001985314 


-0.1402132305 


0.2697081277 


MC 



TABLE II: Same as in Table UJ 
for the last 3 parameters p§, p$, pj of the fit of HNC data. 



p5 


p6 


p7 


Method 


-5.062500000 


-2.0959000000 


0.0676000000 


Ortner 


0.5868725967 


-2.1982700902 


2.7828599024 


HNC - DH 


0.5093239388 


-1.9531860886 


2.5039620685 


HNC 



2. Fits 

We used the functional form of Ortner asymptotic expression 12.11 to fit the HNC data for 
j3u/T in the interval < T < 1. We are left with a eight parameters fit (i.e. the pi 
for % = 0, ... 7) or a seven parameters fit, if p is fixed to its Debye value p = —y/3/2. 
The values found for the pi are given in the Tables H and [Til For the eight parameters fit 
the maximum deviation of the fit from the HNC data is 7.3 10~ 7 with a mean deviation 
of 1.9 10~ 7 , while for the seven parameters fit these deviations are 1.3 10 -6 and 3.3 10 -7 
respectively. Some comments are in order. 

• Firstly, for T < 0.1 the estimations of f3u/T in the framework of HNC, Cohen et al. 
and Ortner theories all coincide with an absolute precision of the order of 1.10 -4 , as 
apparent in table IIHI These conclusions are also true for DH approximation. 

• The agreement between HNC energies and that predicted by Cohen et al. expression 
(cf "Thl" in figure [3] and table HUJ) differ by less than 2.1Cr 3 in the range < V < 0.3. 



Note that the apparent discrepancies between the pi of the fit of HNC and the "exact" 
coefficients of Cohen expansion do not spoil the excellent agreement between the two 
approaches. 

• The agreement between HNC energies and that predicted by Ortner et al. expression 
(cf "Th2" in figure [3] and table [HID differ by less than 2.10" 3 in the range < T < 0.2. 

From these remarks we conclude that HNC is, as expected, exact in the low coupling regime 
at least up to T = 0.3. Moreover DH theory cannot be trusted for T > 0.1, Cohen et al. 
expression can be used confidently as it stands for T < 0.3 and, unexpectedly, the additional 
orders in the asymptotic expression obtained by Ortner do not improve, unfortunately, on 
Cohen results. We suggest to reexamine the details of the calculations of reference 4|]. The 
functional Dependance in T of equation (12. ip is probably correct but misprints in one of the 
Pi for either i = 5,6 or i = 7 are likely. 

C. MC theoretical background 

MC simulations are not well adapted to the low coupling regime for two reasons. First, 
since the configurational energies are small, the convergence of the MC process is slow. 
Secondly, in the case of the OCP considered here, the Debye length Xd = 1/VST diverges 
as T — > and thus becomes larger than the (finite) size of the simulation box, with entails 
severe finite size effects. To use the MC method for obtaining very precise results for the 
OCP in the range of < T < 1 is therefore a real challenge. Some comments on our 
methodology seem to us worthwhile. 

Our simulations were performed in the canonical ensemble within hyperspherical bound- 
ary conditions. The particles are thus confined on the surface of a 4D sphere S 3 of radius 
R and the plasma pair potential between ions is simply the Coulombic interaction in this 
geometry. The latter has a simple analytical expression which allows high precision com- 
putations in contrast with the usual technique of Ewald summations where the potential is 
poorly determined at short distances. The theoretical background of this method has been 



already described in details in previous works [lOj, 111] and will not be rediscussed here. We 
only extract from these previous theoretical considerations the following point. It turns out 
that DH equation (i.e. Helmoltz equation) can be solved analytically in S 3 which yields the 



exact finite size dependence of the excess internal energy in this approximation and therefore 
in the low coupling limit. One finds that at the leading order 



UN 



(r) - Uqo (r) ~ N~ 2/s for r -> and N -> oo . (2.4) 



Of course this behavior in only asymptotic and sub-leading terms in [iV -2 / 3 ] 2 , [iV -2//3 ] 3 
must be taken into account if N is not large enough. For couplings T > 3 we shown in paper 
I that we rather have un (T) — Uoo (T) ~ N~ x . This remark yields the correct procedure : 
for a given parameter V perform MC simulations for different number of particles N and 
take advantage of the scaling relation 12.41 to obtain the thermodynamic limit (T). The 
estimation of the statistical errors on the un (T) and the extrapolated thermodynamic limit 



Moo (r) is also described in details in I. However, by contrast with refs |10l . Ill] devoted to 
the strong correlation regime (1 < T < 190), present work only the small couplings are 
considered. In order to test the validity of HNC, notably in the range 0.3 < T < 1 with 
an error of ~ 1.10 -4 we were led to perform huge Markov chains and consider very large 
systems up to iV = 51200 particles in order to reach the scaling region where [2741 applies. 
Since HNC and Cohen asymptotic forms for u differ by less than ~ 1. 10~ 4 in the range 
< T < 0.3 we can claim (as will be discussed in details below) an overall maximum error 
of ~ 1. 10~ 4 for the dimensionless f3u/T in the whole interval < T < 1. 

Some additional simulations in the transition region to high correlation regime 1 < T < 10 
were also performed to make contact with our previous results. 



III. MC DATA ANALYSIS AND FITS 



A. Data analysis 

We adopted the same procedure as the one described in reference I. The MC simula- 
tions were performed using the standard Metropolis algorithm to build Markov chains in 
the canonical ensemble. In the small T regime, < V < 1, where finite size effects are 
tremendously important, we considered much larger systems than before. In order to get 
the thermodynamic limit (TL) of the excess internal energy for each value of T, we performed 
simulations for samples of iV = 100, 200, 400, 800, 1600, 3200, 6400, 12800, 25600, and 51200 
particles. The cumulated reduced excess energy (CREE) 0U(T, N)/T at coupling T and 



TABLE III: Minus the dimensionless energy (3un(T)/T of the OCP as a function of T for MC (with 
error bars), HNC, Cohen, and Ortner approximations. 



r 


MC 


HNC 


Cohen 


Ortner 


0.1 


0.25117(34) 


0.25688548 


0.25677226 


0.25699174 


0.2 


0.34111(17) 


0.34238929 


0.34127338 


0.34436859 


0.3 


0.397693(64) 


0.39837711 


0.39608173 


0.40761777 


0.4 


0.439323(53) 


0.43968253 


0.44115547 


0.46432208 


0.5 


0.472172(42) 


0.47208481 


0.49302326 


0.521520956 


0.6 


0.498715(21) 


0.49850618 


0.57144385 


0.58565711 


0.7 


0.521064(20) 


0.52064202 


0.70120487 


0.67244493 


0.8 


0.540173(15) 


0.53956586 


0.91276540 


0.81996338 


0.9 


0.556823(30) 


0.55600050 


1.2425017 


1.1053739 


1.0 


0.571403(24) 


0.57045534 


1.7327877 


1.6651877 



number of particles N, was computed as the cumulated mean over M successive configura- 
tions "i" of the Markov chains as 

, M 



(3u N AM) 



E 



M ^ NT 

i=i 



1 < M < n nconf ) , 



(3.1) 



We generated MC chains of n ncon f = 4.10 9 configurations after thermal equilibration, for 
all systems up to N = 25600 particles. The reason was to to reach a stable plateau for 
the CREE and to reduce statistical errors. These two points will be illustrated further. 
For iV = 25600 such long chains result in the mixing of 5 independent chains, each one 
corresponding to half a month of CPU time. Thus the N = 25600 value of the excess 
energy represent a 2 months and a half calculation. For N = 12800 the total duration was 2 
monthes, with two independent chains. For comparison a N = 800 simulation is performed 
in 2 days in a unique chain. One day is enough for aJV = 400 simulation. These calculations 
have been performed simultaneously on the CEA Opteron clusters, local PC and the CRI 
cluster of Orsay, using one processor by job. 

In order to compute MC statistical errors on f3u^{T)/T each total run was divided into 
ub blocks and the error bar was obtained by a standard block analysis 17| . Each block 
involved a large number n^ 1 ^ of successive MC configurations and was supposed to be 



TABLE IV: Minus the MC energy (3un(T)/T of the OCP as a function of T and the number of 
particles N. The number in bracket which corresponds to one standard deviation a is the accuracy 
of the last digits. 



r 


N = 1600 


N = 3200 


N = 6400 


N = 12800 


N = 25600 


N = 51200 


0.1 


0.20942(7) 


0.21343(8) 


0.21984(8) 


0.22728(7) 


0.23436(7) 


0.24018(18) 


0.2 


0.32066(7) 


0.32477(6) 


0.32865(4) 


0.33223(4) 


0.33502(4) 


0.337244(74) 


0.3 


0.385447(36) 


0.388389(23) 


0.391016(23) 


0.393158(27) 


0.394704(25) 


0.395825(99) 


0.4 


0.430965(23) 


0.433233(23) 


0.435009(19) 


0.436452(19) 


0.437507(18) 


0.438230(63) 


0.5 


0.465821(16) 


0.467568(17) 


0.468939(17) 


0.470015(17) 


0.470754(15) 


0.471263(52) 


0.6 


0.493812(13) 


0.495220(13) 


0.496352(16) 


0.497158(13) 


0.497714(13) 


0.498072(37) 


0.7 


0.517109(13) 


0.518254(12) 


0.519178(12) 


0.519806(11) 


0.520259(13) 


0.520600(32) 


0.8 


0.536909(8) 


0.537854(7) 


0.538606(9) 


0.539140(11) 


0.539513(11) 


0.539745(25) 


0.9 


0.554034(10) 


0.554810(12) 


0.555458(8) 


0.555886(11) 


0.556232(10) 


0.556458(40) 


1.0 


0.569012(15) 


0.569714(9) 


0.570281(8) 


0.570669(10) 


0.570930(9) 


0.571119(24) 



statistically independent of the others. For each calculation we checked that the variance 
was independent of the size of the blocks for sufficiently large values of n c ^ n ^ . Results are so 
stable that we shall no more discuss this point in this paper. The need of large simulations 
with iV = 51200 particles appeared with the difficulty to reach the thermodynamic limit 
and to obtain the wanted precision for the T values that we considered. But, due to huge 
demand in CPU time of these simulations (one month for 10000 configurations) only short 
chains were considered, however long enough to reach the stable plateau of the CREE and to 
improve the TL research (see below). Our data for /?Wjv(r)/r are reported in table HVl where 
the number in bracket correspond to one standard deviation a and represent the accuracy 
of the last digits and only the results for iV > 1600 are given. 

B. Connection with former simulations for 1 < T < 10 

Before we present our new results for < T < 1, we shall study the connection with 
the results obtained in paper I, calculated with the same MC code, but another range 
of r values, i.e. T > 1. The only difference between the two calculations, calculated in 



TABLE V: Comparison with previous results of the MC energy /3un(T)/T of the OCP in space 5 3 
in function of the number of particles N for T = 5 and T = 10. The first row, case "a" , corresponds 
to present study and second row, case "b", to table 1 of The only difference between the 



two calculations is the number of configurations, typically n con f — 800 10 6 MC configurations after 
equilibration in case " a" , and n con f = 5.10 9 in case "b". The number in bracket which corresponds 
to one standard deviation a is the accuracy of the last digits. With two standard deviations the 
agreement is fulfilled. 



r 


N = 400 


N = 800 


N = 1600 


N = 3200 


N = 6400 


case 


5 




.7510930(37) 


.7511501(31) 


.7512037(35) 


.7512332(21) 


a 


5 


.7510201(89) 


.7511042(126) 


.7511513(135) 


.7511775(85) 




b 


10 




.7998396(26) 


.7998148(30) 


.7998098(28) 


.7998043(15) 


a 


10 


.7998865(53) 


.7998414(43) 


.7998149(51) 


.7998131(55) 




b 



double precision on 64 bytes work stations, is thus the maximum number of configurations, 
typically n conf = 800 10 6 MC configurations -after equilibration in previous case (case "a"), 
and n con f = 5.10 9 in this paper (case "b" ). We have performed comparisons for T = 1, 2, 3, 4, 5 
and T = 10. The choices retained in I were, at that time, the maximum reasonable conditions 
for the simulations. 

Finite size effects decrease with increasing value of T. Details of CREE's for T = 5 
and r = 10 are reported in Table [V] (N dependence of MC energy Pu^(T)/T in both 
calculations, present and I). For T = 10 the results are in good agreement. But for T = 5 
slight discrepancies observed at iV = 1600 and N = 3200 between the two calculations are a 
bit worrying. Indeed in these cases the error bars intervals do not overlap. The main reason 
is that, for the lowest T results of ref. I. the plateau of the CREE was in fact not reached. 
This feature is illustrated by figure [5] where the CREE's for T = 2 are displayed. The figure 
illustrates the lack of configurations in the simulations of ref. I for the CREE /3U/T versus the 
number of configurations, displayed for different number of particles. From top to bottom 
N = 800, 1600, 3200, 6400. The blue arrow points on the maximum number of configurations 
considered in I. When compared to our new calculations, clearly the Markov chain was not 
long enough to reach a plateau and such a drift of the CREE was probably underestimated in 
our previous calculations. The large variation with N of the CREE with N ( solid red line ) 
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FIG. 5: Solid lines: cumulated reduced excess energy f3u]\f(T) /T versus the number of configurations 
for r = 2. From bottom to top N = 800,1600,3200,6400. Symbols: block averages. The blue 
arrow points to the maximum number of configurations n con f = 8. 10 s considered in ref. I. 

gives an idea of the amplitude of finite size effects. The simulation for the case N = 6400, not 
included in paper I, has been added to improve the TL extrapolation. Only the 4. 10 9 first 
configurations are plotted for visibility, but clearly each CREE value reaches its equilibrium 
value for a fixed iV value. Figure [6] illustrates how previous conclusions for the case T = 2 
are emphasized in the case T = 1. It follows from the above remarks that a re- analysis of 
the TL of the energy of the OCP is necessary for T = 1, 2, ... 10. We recall the conclusions 
of I according to which the scaling law 12.41 is valid only for low T and that for T > 3 — 4 the 
thermodynamic limit is reached more quickly with a scaling law 



u N (r) - Moo (r) ~ N- 1 for r > 3 - 4 and iV -> oo . 



(3.2) 



Moreover the scaling limits l2~4l and I3~2l are satisfied, depending on the value of T, for very 
large, and sometimes prohibitive large, numbers of particles N. The ideal case would be 
a linear fit passing through all the points within the error bars. This situation was indeed 
observed by including simulations at iV = 51200 particles and for not too low values of T. In 
other situations we had to content ourselves with quadratic fits including the next leading 
order term (i.e. either 0(1/ N 2 ) or (9(1 /TV 4 / 3 ) according to the value of T). In Table PVTl are 
resumed the comparisons for the TL of the energy between present and results of paper I 
for T = 1, 2, 3, 4, 5, 10. The type of the extrapolation scheme is specified in the column "fit", 
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FIG. 6: Solid lines: cumulated reduced excess energy f3uj^(Y) /Y versus the number of configurations 
for T = 1. From top to bottom N = 3200, 6400, 12800, 25600, 51200. Symbols: block averages. 

together with the interval of iV values considered for the fit. Precision are also reported. 
For T < 5 it is clear that results are slightly shifted between calculations. As expected for 
r = 10 present and previous results are similar; higher values of Y should not cause any 
trouble. 



C. Thermodynamic limit extrapolation scheme 

The aim of our simulations was to compute the TL of the energy /3un =qo (Y) with a 
high degree of accuracy by taking into account finite size effects which are of overwhelming 
importance for T < 1. The need of simulations up to iV = 51200 and involving no less than 
N = 800, or even N = 1600 particles for the smallest values of Y, appeared crucial to reach 
the scaling law 12.41 It appears that, for this range of N, MC data can be fitted with the 
quadratic fits 

(3.3) 

For most values of Y it proved possible to explicitely check the asymptotic form linear 
in TV -2 / 3 (i.e. ai = in equation 13. 3ft by keeping only the 3 largest systems, i.e. N = 
12800, 25600 and N = 51200. Recall that in paper I the largest considered systems were 
made of N = 3200 particles. An exhaustive discussion follows in next section. 



1 1 

(3u N {Y) = (3u N=OQ {Y) + axj^ + a 2 



TABLE VI: Thermodynamic limit of the energy of the OCP versus T for T > 1 of , case "a", com- 
pared to previous calculations, case "b". The difference between two calculations is the maximum 
number of partparticles (no more than 3200 in case a) and the total number of configurations. The 
type of extrapolation scheme is specified in the column "fit". For instance quad(3200 — 51200) 
means that a quadratic regression involving the data from N = 3200, 6400, 12800, 51200 has been 
used. The variable entering the fit is specified in the column Variable, p is the precision of the fit. 
The number in bracket which corresponds to one standard deviation a is the accuracy of the last 
digits. 



r 




Fit 


p* 10 5 


/3iioo/r 


Fit 


p* 10 5 


Variable 




a 


a 


a 


b 


b 


b 




i. 


-0.571387(24) 


lin(12800-51200) 


4.2 


-0.571098(39) 


cub(100-3200) 


6.9 


AT-2/3 


i. 


-0.571403(22) 


quad(3200-51200) 


3.8 








AT-2/3 


2. 


-0.6598934(68) 


quad(800-6400) 


7.0 


-0.659983(23) 


quad(200-3200) 


3.5 


AT-2/3 


3. 


-0.7042987(54) 


quad(3200-6400 


0.8 


-0.704348(19) 


quad (200-3200) 


2.7 


AT-2/3 


4. 


-0.7319760(46) 


lin(800-6400) 


0.6 


-0.731916(12) 


quad (200-3200) 


1.7 


N- 1 


5. 


-0.7512608(22) 


lin(800-6400) 


0.3 


-0.7512126(98) 


quad (200-3200) 


1.3 


AT-i 


10. 


-0.7997991(16) 


lin(800-6400) 


0.2 


-0.7997974(45) 


lin(400-3200) 


0.56 


at- 1 



D. Results for < T < 1 

We present and discuss in details the ten values T = 0.1, 0.2, . . . , 1 considered in our numer- 
ical experiments. Figures 171 IB| 191 and [TOl illustrate the CREE /Jujv(r)/T versus the number 
of configurations for several caracteristic values of T( T = 0.1, 0.2, 0.4, and 0.7 respectively) 
typical of the different plasma regimes in the interval (0,1). Our previous comments on 
figures [6] and [5] (for F = 1, 2 respectively) are still valid in these cases. We stress once again 
the need of large systems together with the need of enough configurations to reach a stable 
plateau after thermal equilibration. 

All generated configurations, n con f = 6.10 9 , are displayed in figure [7] (r = 0.1) while 
a zoom of only the first 2. 10 9 configurations is displayed in figure [8] (r = 0.2), which 
exemplifies the plateau reached by the CREE for A" = 51200. On the last two figures [9] and 
[101 respectively for T = 0.4 and T = 0.7 and n con f = 4. 10 9 , we see the good convergence with 



TABLE VII: Thermodynamic limit of the energy of the OCP versus V. The number in bracket 
which corresponds to one standard deviation a is the accuracy of the last digits. The type of 
extrapolation scheme is specified in the column "fit". The variable entering the fit is A -2 / 3 . 
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/3uoo/r 


Fit 


0.1 


-0.25117(34) 


quad(6400-51200) 


0.2 


-0.34111(17) 


quad(6400-51200) 


0.3 


-0.397693(64) 


quad(3200-51200) 


0.4 


-0.439323(53) 


lin(12800-51200) 


0.4 


-0.439528(50) 


quad(3200-51200) 


0.5 


-0.472028(45) 


lin(12800-51200) 


0.5 


-0.472172(42) 


quad(3200-51200) 


0.6 


-0.498663(38) 


lin(12800-51200) 


0.6 


-0.498715(21) 


quad(1600-51200) 


0.7 


-0.521063(32) 


lin(12800-51200) 


0.7 


-0.521064(20) 


quad(1600-51200) 


0.8 


-0.540146(28) 


lin(12800-51200) 


0.8 


-0.540173(15) 


quad(1600-51200) 


0.9 


-0.556823(30) 


lin(12800-51200) 


0.9 


-0.556801(25) 


quad(3200-51200) 


1.0 


-0.571387(24) 


lin(12800-51200) 


1.0 


-0.571403(22) 


quad(3200-51200) 



A" as the interval width between CREE values decreases from top to bottom. By contrast 
the low T runs do not exhibit this regular decrease. Of course beyond visual impressions 
only the possibility and precision of the fitting process of the MC CREE results will give a 
firm answer on the quality of the TL calculation for each T value. 

Table HVl resumes present work MC calculations of the MC energy /?ttjv(r)/r of the OCP 
as a function of T for A" = 1600 to N = 51200. The number in bracket which corresponds 
to one standard deviation a is the accuracy of the last digits. Results corresponding to 
A" < 1600, not included in the fits, are not reported. The thermodynamic limit values of 
the energy versus V are reported in Table IVHI The type of extrapolation schemes retained 



in the fits, i.e. linear or quadratic (cf equation I3.3[) . are specified. In figures [TT j [T2l [TBI and 
[14] we display the quadratic fit of /?Mjv(r)/T (solid black line) and the linear fits for the 3 
largest numbers of particles considered, when available (red dashed line) for T = 1, 0.7, 0.4, 
and T = 0.1 respectively. The error bars on the value of the TL of the energy Puoo(T)/r 
reported in table IVHI are the error bars of the linear (or quadratic) regression. 

We discuss now the results from high to low T's. Figure [11] illustrates the high quality 
of the fits obtained in the case T = 1.0. Indeed the extrapolated TL of /3ujv(r)/r coincide 
for both the linear and the quadratic fits (the latter involving more states with low number 
of particles and the former only the 3 largest systems) with a nice overlap of the error bars. 
Results are similar down to V = 0.7, as illustrated in Figure [121 for T = 0.7. 

In the range 0.5 < T < 0.7 the precision of the fits is good but the linear and the quadratic 
fit extrapolations do not give exactly the same TL values, however the error bars do overlap. 
Figure [13] corresponding to the case T = 0.4, illustrates the smallest T at which a linear fit 
is possible with the 3 higher values of N. The linear and the quadratic fit extrapolations 
giving the TL values would coincide within the error bars if the latter were defined to be 
two standard deviations rather than only one according to our choice. 

For T smaller or equal to 0.3 it was impossible to reach an asymptotic form of /3un(T)/T, 
linear in the variable N~ 2 ^ 3 , and only a quadratic polynomial fit was possible (cf table 
IVIIj) . For that reason it is legitimate to consider the error bars on the extrapolated value 
/3iioo(r)/r as overoptimistic in this range of T, see figure [14] for an illustration in the case 
T = 0.1. Simulations involving larger numbers of particles would be necessary but are out 
of our reach. 

For all the states with a T > 0.4 the TL limit Woo(r) can thus be obtained with a high 
precision p ~ 10~ 5 , after a careful study of finite size effects on the MC energies un{F). For 
smaller values of T, for instance T = 0.1, samples of more than N ~ 200000 particles should 
be used to reach the leading order of the asymptotic regime 12.41 However such an effort 
would be useless since HNC and Cohen approximations are then "exact" within the wanted 
precision on u. The Woo(r) are perfectly well fitted in the range 0.4 < V < 1 by the Cohen's 
functional form, given by equation 12. lal involving the five parameters Pi (i — 0, . . . , 4) given 
m table H 



IV. CONCLUSION 



In this conclusion we compare at first the Cohen-Ortner low T expansions, HNC and MC 
data. Figure [3] shows without ambiguity the good agreement between HNC and MC results 
in the range < T < 1 and the large departure of both results with analytical expansion 
ones, DH (for T > 0.05), Thl (for T > 0.3) and Th2 (for T > 0.2). Note however that 
the scale of the figure is not large enough to discriminate between HNC and MC results, 
notably because the errors bars on MC results are smaller than the size of the symbols. A 
more enlightening illustration is that of figure [2] which gives the ratio of the MC and HNC 
energies. The disagreement for T < 0.3 results from a bad evaluation of the TL of due 
to huge finite size effects spoiling the MC data, while the disagreement for V > 0.6 simply 
reflects the failure of the HNC approximation at high V . A nearly perfect agreement between 
MC and HNC results, compatible with one standard deviation is observed only at V = 0.5; 
with two standard deviations the HNC results are within the error bars of the MC data in 
the interval 0.4 < T < 0.6. By passing our new HNC calculations for some values of V in 
the range (1, 10) are plotted in figure H] were MC data were also included for comparison. 

It is the place to resume our analysis. We found that, for a wanted precision of p = 10~ 5 
on the energy : 

• < T < 0.05 is the range of validity of Debye-Hiickel theory. 

• < T < 0.3 is the range of validity of Cohen low T expansion 12.11 

• Ortner's additional terms do not improve the results. 

• < T < 0.5 is the range of validity of HNC. The data are perfectly represented by the 
eight parameters fit of tables fl] and [Til 

• We were able to extract the thermodynamic limit of the OCP energy from our MC 
simulations with a precision not smaller than p = 10 -5 in the range 0.4 < T < 1. Our 
data are well fitted by the five parameters fit of table [TJ 
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FIG. 8: Solid lines: cumulated reduced excess energy Pun(T) /T versus the number of configurations 
for r = 0.2. From top to bottom N = 400, 800, 1600, 3200, 6400, 12800, 25600, 51200. Symbols: 
block averages. 
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FIG. 9: Solid lines: cumulated reduced excess energy f3uj^(Y) /Y versus the number of configurations 
for r = 0.4. From top to bottom N = 400, 800, 1600, 3200, 6400, 12800, 25600, 51200. Symbols: 
block averages. 
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FIG. 10: Solid lines: cumulated reduced excess energy /3un(Y)/Y versus the number of config- 
urations for r = 0.7. From top to bottom N = 400,800,1600,3200,6400,12800,25600,51200. 
Symbols: block averages. 




FIG. 11: Solid lines: cumulated reduced excess energy /3uAr(T)/r versus I/N 2 / 3 for T = 1. From 
left to right N = 00, 51200, 25600, 12800, 6400, 3200. The error bars correspond to one standard 
deviation a. Solid black line : quadratic polynomial regression of MC data. Dashed red line : 
linear regression of the 3 larger systems MC data. 




FIG. 12: Same legend than figure [TT] but for T 
00, 51200, 25600, 12800, 6400, 3200, 1600. 



= 0.7. From left to right N = 




FIG. 13: Same legend than figure [IT] but for T = 0.4. From left to right N = 
oo, 51200, 25600, 12800, 6400, 3200. 
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FIG. 14: Solid lines: cumulated reduced excess energy f3uN(T)/T versus 1/iV 2 / 3 for T = 0.1. From 
left to right N = oo, 51200, 25600, 12800, 6400, 3200. The error bars correspond to one standard 
deviation a. Solid black line : quadratic polynomial regression of MC data. 



